#include "../lj.h"
#include "../ini/velocity.h"
#include "../sys/rand.h"
#include "../sys/wtime.h"
#include <time.h>
#include "../io/xyz.h"
#include "../io/gro.h"
#include "../step/verlet.h"
#include <math.h>
#include "../constants.h"
#include "../ini/crystal.h"
#include "../dynamics.h"

#define ERROR program_result=-1; goto END


int main( int argc, char *argv[] )
{
	
	int              program_result = 0;
	mds_lj_data      lj;
	mds_frame        frame    = MDS_FRAME_INIT;
	mds_vector_set   position = NULL;
	mds_vector_set   velocity = NULL;
	size_t           np       = 0;
	Real             box_size = 0;
	const Real       epsilon  = 0.995708; /*-- for argon */
	const Real       sigma    = 0.3405;   /*-- for argon */
	mds_step         step;
	mds_verlet_data  verlet_data;
	Real             dt = 0.0;
	Real             t_min=0;
	FILE            *fp   = NULL;
	Real             temperature = 100;
	Real             K0          = 0;
	mds_uvector      cr;
	mds_dynamics     dynamics = { mds_dynamics_full, NULL };
	register size_t  i;
	Real mass_min = 1e9;
	
	mds_srand( time(NULL) );
	
	if( argc <= 2 )
	{
		fprintf( stderr, "usage: %s num_particles box_size\n", argv[0]);
		ERROR;
	}
	
	np       = strtol( argv[1], NULL, 10 );
	box_size = strtod( argv[2], NULL);
	
	frame.box.x = box_size;
	frame.box.y = box_size;
	frame.box.z = box_size;
	fprintf( stderr, "md> query num_particles=%u\n", (unsigned)np );
	{
		np = mds_crystal_close_packed_count( &cr, np, & frame.box);
	}
	fprintf( stderr, "md> got   num_particles=%u\n", (unsigned)np );
	
	
	/** creating frame */
	if( mds_frame_load( &frame, np, MDS_NUM_SETS_FOR_VERLET) < 0 )
	{
		mds_perror( stderr, "can't load frame(%u,%u)\n", (unsigned) np, MDS_NUM_SETS_FOR_VERLET); 
		ERROR;
	}
	
	/** finalize frame */
	position = frame.sets[ MDS_POSITION ];
	velocity = frame.sets[ MDS_VELOCITY ];
	frame.pbc_for_f = mds_pbc_xyz;
	frame.pbc_for_r = mds_pbc_xyz;
	
#if 0
	 frame.pbc_for_f = mds_pbc_none;
	 frame.pbc_for_r = mds_pbc_none;
#endif
	
	/* Real ix[frame->num_particles]; */
	
	
	
	/** initialize masses */
	fprintf( stderr, "md> initialize masses\n" );
	for( i=0; i < np; ++i )
	{
		const Real mass = 10 + 10 * mds_alea();
		position[i].info.inv_mass = 1.0 / mass;
		velocity[i].info.mass     = mass;
		if( mass < mass_min ) mass_min  = mass;
	}
	
	fprintf( stderr, "md> min mass=%g \n", mass_min );
	
	/** initialize position */
	fprintf( stderr, "md> initialize positions\n");
	/* mds_frame_initialize_positions( &frame ); */
	mds_crystal_build_fcc(position, np, &cr, & frame.box, 0);
	
	/** initialize velocities */
	fprintf( stderr, "md> initialize velocities\n");
	mds_frame_initialize_velocities( &frame, temperature);
	K0 = frame.K;
	fprintf( stderr, "md> temperature= %g K, K0 = %g\n", temperature, K0 ); 
	
	/** initialize potential parameters */
	mds_lj_init( &lj, epsilon, sigma );
	
	
	/** compute approx min LJ harmonic period */
	t_min = sigma * sqrt( mass_min / epsilon );
	fprintf( stderr, "t_min= %g ps\n", t_min );
	t_min  *= 0.0005;
	dt     = 1;
	while( dt > t_min ) dt *= 0.1;
	fprintf( stderr, "dt  s= %g ps\n", dt );
	fprintf( stderr, "T^* = %g\n", MDS_R * temperature / epsilon);
	
	
	
	/** Initialize Accelerations */
	mds_dynamics_call( &dynamics, &frame, mds_lj, & lj );
	
	/** load a step */
	mds_nothing_step_load( &step, NULL);
	mds_step_init_dt( &step, &frame, dt);
	
	/** load verlet step */
	mds_verlet_step_load( &step, &verlet_data );
	mds_step_init_dt( &step,&frame, dt );
	
	fprintf( stderr, "K= %g kJ/mol\n", frame.K );
	fprintf( stderr, "V= %g kJ/mol\n", frame.V );
	fprintf( stderr, "W= %g kJ/mol\n", frame.W );
	
	
	
	mds_xyz_save("md.xyz", position, np, "first frame", 0);
	fp = fopen( "info.txt", "wt" );
	if( fp )
	{
		fprintf( fp, "0 %g %g %g %g\n", frame.V, frame.K, frame.K + frame.V, frame.W );
		fclose(fp);
	}
		
	{
		Real     t     = 0;
		unsigned istep = 0;
		Real     old_t = 0;
		
		while( t < 20.0 )
		{
			mds_step_predict( &step, &frame );
			mds_dynamics_call( &dynamics, &frame, mds_lj, &lj );
			mds_step_correct( &step, &frame);
			
			if( t < 10 )
			{
				frame.K = mds_velocity_rescale( velocity, np, sqrt(K0/frame.K) );
			}
			
			++istep;
			t = istep * dt;
			
			
			if( (t-old_t) >= 0.05 )
			{
				fprintf( stderr, "t= %8g | V = %8.3e | K = %8.3e | W = %8g | E = %8.3e        \r", t, frame.V, frame.K, frame.W, frame.K + frame.V );
				mds_xyz_save("md.xyz", position, np, "other frame", 1);
				fp = fopen( "info.txt", "at" );
				if( fp )
				{
					fprintf(fp, "%g %g %g %g %g\n", t, frame.V, frame.K, frame.K + frame.V, frame.W );
					fclose(fp);
				}
				old_t = t;
			}
			
		}
		fprintf( stderr,"\n\n");
		
	}
	
	
	
	
END:
	MDS_MEMORY();
	mds_frame_quit( &frame );
	MDS_MEMORY();
	return program_result;
}
